
##############################################################
## Replication code
## Alsaadi, "Unconditional Loyalty: The Survival of Minority Autocracies"
##############################################################

## Packages ###########################################################################################################
library(tidyverse)
library(survival)
library(foreign)
library(survminer)
library(cobalt)
library(cowplot)
library(marginaleffects)
library(lmtest)
library(stargazer)
library(sandwich)
library(MatchIt)
library(optmatch)
library(mitools)
library(Amelia)
library(Synth)
library(gsynth)
library(patchwork)
library(ggpubr)





#############################################################################################################


# Load data 
minority_TSC=readRDS("minorityreg.RDS")

names(minority_TSC)
###########################
# Descriptive statistics in the introduction 
###########################
#On average, minority regimes with excluded majority lasted more than twice as long as other autocracies
Average_tenure<-minority_TSC %>% 
  group_by(cus_caseid) %>% 
  dplyr::select(cus_t_surv, minority) %>% 
  summarise_all(mean)

Average_tenure %>% group_by(minority) %>% 
  summarise(tenure= mean(cus_t_surv)) #### (34 for minority regimes and 15 for non-minority regimes)


###########################
# Figure 1. Survival Curves: Autocratic Regimes (1900–2015)
###########################
surv_min <-survfit(Surv(cus_t,cus_t+1,cus_fail) ~ minority, data=minority_TSC)
surv_allmin <-survfit(Surv(cus_t,cus_t+1,cus_fail) ~ allminority, data=minority_TSC)
surv_fracmin<- survfit(Surv(cus_t,cus_t+1,cus_fail) ~ frac_minority, data=minority_TSC)


surv1=ggsurvplot(surv_allmin , data=minority_TSC,
                 palette = c("#2F2F2F", "#979797"),
                 risk.table=FALSE,
                 conf.int = TRUE,
                 xlim=c(0,110),
                 break.time.by =10, 
                 ggtheme = theme_minimal(base_size = 14) +   
                   theme(
                     axis.title.x = element_text(size = 19),    
                     axis.title.y = element_text(size = 19),    
                     legend.title = element_text(size = 19),    
                     legend.text = element_text(size = 19),
                     plot.margin = margin(20, 20, 20, 20) 
                   ),
                 risk.table.y.text=FALSE,
                 censor.size =1,
                 legend.title="a",
                 legend.labs=c("Other regimes","All minority regimes"))


surv2= ggsurvplot(surv_min , data=minority_TSC,
                  palette = c("#2F2F2F", "#979797"),
                  risk.table=FALSE,
                  conf.int = TRUE,
                  xlim=c(0,110),
                  break.time.by =10, 
                  ggtheme = theme_minimal(base_size = 16) +   
                    theme(
                      axis.title.x = element_text(size = 19),    
                      axis.title.y = element_text(size = 19),    
                      legend.title = element_text(size = 19),    
                      legend.text = element_text(size = 19),
                      plot.margin = margin(20, 20, 20, 20) 
                    ),
                  risk.table.y.text=FALSE,
                  censor.size =1,
                  legend.title="b",
                  legend.labs=c("Other regimes","Minority excluding majority"))



surv3= ggsurvplot(surv_fracmin, data=minority_TSC,
                  palette = c("#161616", "#979797"),
                  risk.table=FALSE,
                  conf.int = TRUE, 
                  xlim=c(0,100),
                  break.time.by =10, 
                  ggtheme = theme_minimal(base_size = 14) +   
                    theme(
                      axis.title.x = element_text(size = 19),    
                      axis.title.y = element_text(size = 19),    
                      legend.title = element_text(size = 19),    
                      legend.text = element_text(size = 19),
                      plot.margin = margin(20, 20, 20, 20) 
                    ),
                  risk.table.y.text=FALSE,
                  censor.size =1,
                  legend.title="c",
                  legend.labs=c("Other regimes", "Other minority regimes"))


fig1=ggarrange(surv1$plot,                                                 
          ggarrange(surv2$plot, surv3$plot, ncol = 2),
          nrow = 2)

pdf("fig1.pdf",width=13,height=10)
print(fig1)
dev.off()


# Log rank test
diff.out <- survdiff(Surv(cus_t_surv,cus_fail) ~ allminority, data =minority_TSC) # p= 0.02


# 25-year survival rates
km_min <- survfit(Surv(cus_t,cus_t+1,cus_fail) ~ 1, data = minority_TSC, subset=minority==1)
survest_min <- stepfun(km_min$time, c(1, km_min$surv))
survest_min(25) # 58% for minority regimes 

km_min2 <- survfit(Surv(cus_t,cus_t+1,cus_fail) ~ 1, data = minority_TSC, subset = minority == 0)
survest_min2 <- stepfun(km_min2$time, c(1, km_min2$surv))
survest_min2(25) # 25% for non-minority regimes





###########################
# Regression models
###########################

# Table 2: Authoritarian Regimes Breakdown (1900-2015)
modf1= glm(cus_fail~ allminority, data=minority_TSC, family="binomial")
modf2= glm(cus_fail~ allminority+ log_e_gdppc+ gdp_growth+ log_e_pop+ log_oil_income_pc+ partB+ mon+ persB+mil+
              foreignsupport+ factor(Country)+ factor(Year), data=minority_TSC, family="binomial")

modf3= glm(cus_fail~ minority, data=minority_TSC, family="binomial")
modf4= glm(cus_fail~ minority+ log_e_gdppc+ gdp_growth+ log_e_pop+ log_oil_income_pc+ partB+ mon+ persB+mil+
              foreignsupport+ factor(Country)+ factor(Year), data=minority_TSC, family="binomial")

modf5= glm(cus_fail~ frac_minority, data=minority_TSC, family="binomial")
modf6= glm(cus_fail~ frac_minority+ log_e_gdppc+ gdp_growth+log_e_pop+ log_oil_income_pc+ partB+ mon+ persB+mil+
              foreignsupport+ factor(Country)+factor(Year),
           data=minority_TSC, family="binomial")

#Cluster Standard Errors by country
modf1=coeftest(modf1, vcovCL, cluster= minority_TSC$Country)
modf2=coeftest(modf2, vcovCL, cluster= minority_TSC$Country)
modf3=coeftest(modf3, vcovCL, cluster= minority_TSC$Country)
modf4=coeftest(modf4, vcovCL, cluster= minority_TSC$Country)
modf5=coeftest(modf5, vcovCL, cluster= minority_TSC$Country)
modf6=coeftest(modf6, vcovCL, cluster= minority_TSC$Country)


stargazer( modf1,modf2, modf3,modf4, modf5,modf6, omit = c("Country","Year"), type = "text",
           title = "Authoritarian Regime Breakdown",
           covariate.labels = c("All minority regime", "Minority excluding majority","Minority excluding minorities", "Log GDP","GDP growth", "Log population", 
                                "Log oil wealth",
                                "Party", "Monarchy","Personalist","Military", "Foreign support"),
           dep.var.caption = "Dependent Variable: Regime Breakdown",
           no.space = TRUE,  # Remove space between columns
           notes = c("Standard errors clustered by country"),
           add.lines=list(c("Country and Year FE", "No", "Yes", "No", "Yes","No","Yes")),
           model.names = FALSE,
           dep.var.labels.include = FALSE,
           order = c("allminority", "minority", "frac_minority", "log_e_gdppc", "gdp_growth", 
                     "log_oil_income_pc", "log_e_pop", "partB", "mon", "persB","mil",
                      "foreignsupport" ))


modf3= glm(cus_fail~ minority, data=minority_TSC, family="binomial")
p=predictions(modf3, newdata = datagrid(minority = c(0, 1)),vcov = TRUE )
round(p[2,2]*100, digits = 2)
round(p[1,2]*100, digits = 2) # (2.1% vs. 5.6%)





###########################
# Matching
###########################
# Load data 
dat <- readRDS("dat.RDS")

### multiple imputations
imputation_vars_ts <- c("minority", 'log_e_pop', "log_e_gdppc", "prev_partB", "prev_persB",
                        "prev_mon" ,"e_civil_war", 'decolonize',
                        "foreignsupport", 'Year', 'ccode', 'cyear')


m <- 20
set.seed(1237)
d_mi_ts <- amelia(dat[, imputation_vars_ts], ts = 'Year', cs = 'ccode', idvars = 'cyear', m = m, p2s = 1)


additional_vars <- c('Country', 'cus_t', 'cus_caseid', 'cus_t_surv', 'cus_fail')

d_mi_ts$imputations <- lapply(d_mi_ts$imputations, function(x) {
  for (var in additional_vars) {
    x[[var]] <- dat[[var]]
  }
  x[['cus_tplus1']] <- x[['cus_t']] + 1
  x[['cus_caseid']] <- as.character(x[['cus_caseid']])
  return(x)
})



# Combine imputed datasets
dd.mi <- NULL
for (i in 1:d_mi_ts$m) {
  dd <- d_mi_ts$imputations[[i]]
  dd <- dd[!is.na(dd$cus_caseid) & dd$cus_caseid != "", ]  # Keep only observations with content in cus_caseid
  dd$.imp <- i  # Add imputation indicator
  dd.mi <- rbind(dd.mi, dd)
}



# lists to store results
match_results <- list()
mod_fail_results <- list()


for (i in 1:d_mi_ts$m) {
  set.seed(1237 + i)
  dd <- d_mi_ts$imputations[[i]]
  dd <- dd[!is.na(dd$cus_caseid) & dd$cus_caseid != "", ]
  m.out <- matchit(minority ~ log_e_gdppc + log_e_pop + foreignsupport + e_civil_war+
                     decolonize + prev_partB + prev_mon + prev_persB,
                   data = dd, method = "cem", distance = "glm", link = "logit", estimand = "ATT")
  data_match <- match.data(m.out)
  data_match$treat <- data_match$minority
  mod_fail <- glm(cus_fail ~ minority, data = data_match, weights = weights, family = "binomial")
  
  match_results[[i]] <- m.out
  mod_fail_results[[i]] <- mod_fail
}


pooled_mod_fail <- MIcombine(mod_fail_results)
summary(pooled_mod_fail)

# Marginal effects
compute_avg_comparisons <- function(model, data) {
  comparisons <- avg_comparisons(model,
                                 variables = "minority",
                                 vcov = ~subclass,
                                 newdata = subset(data, treat == 1),
                                 wts = "weights")
  return(comparisons)
}


avg_comparisons_list <- lapply(mod_fail_results, compute_avg_comparisons, data = data_match)
ame_list <- do.call(rbind, lapply(avg_comparisons_list, function(x) as.data.frame(x)))
ame_list <- ame_list %>% filter(term == "minority")

pooled_ame <- ame_list %>%
  group_by(term) %>%
  summarize(
    AME = mean(estimate, na.rm = TRUE),
    SE = sqrt(mean(std.error^2, na.rm = TRUE) + var(estimate, na.rm = TRUE) / length(estimate)),
    .groups = "drop"
  )


# Covariate Balance
newNames <- c(
  'log_e_gdppc' = 'GDP per capita (log)',
  'log_e_pop' = 'Population (log)',
  'foreignsupport' = 'foreign/colonial support',
  'e_civil_war'='Civil War',
  'decolonize' = 'Decolonization',
  'prev_partB' = 'Previous party',
  'prev_mon' = ' Previous Monarchy',
  'prev_persB' = ' Previous Personalist'
)

love.plot(m.out, thresholds = c(m = 0.1),
          threshold = 0.1,
          stat = 'mean.diffs',
          limits = c(-1, 1),
          var.names = newNames,
          position = "bottom",
)






###########################################################################################################


## Mechanisms of Authoritarian Survival ###########################################################################################################
# Load NAVCO Data 
navco_autocracy2013 <- readRDS("navco_autocracy.RDS")

# Table 3: Challenges to Authoritarian Regimes (1945-2013)
mod3= glm(success~ minority, data=navco_autocracy2013, family="binomial")
mod3.2= glm(success~ minority+ log_e_gdppc+ log_e_pop+ camp_size_cat+
              regime_support+ camp_support+ camp_confl_intensity+
              prim_meth+log_duration_days+
              factor(location)+ factor(year),
            data=navco_autocracy2013, family="binomial") 

def3= glm(state_defect~ minority, data=navco_autocracy2013, family="binomial") 

def3.2= glm(state_defect~ minority+log_e_gdppc+log_e_pop+ camp_size_cat+
              regime_support+camp_support+
              camp_support+camp_confl_intensity+ prim_meth+log_duration_days+
              factor(location)+factor(year), data=navco_autocracy2013, family="binomial") 

mob3= lm(v2caautmob~ minority, data=navco_autocracy2013) 
mob3.2= lm(v2caautmob~ minority+log_e_gdppc+log_e_pop+ camp_size_cat+ 
             regime_support+camp_support+
             camp_support+camp_confl_intensity+ prim_meth+log_duration_days+
             factor(location)+factor(year)
           , data=navco_autocracy2013) 

#Cluster Standard Errors by country
mod3=coeftest(mod3, vcovHAC, cluster= navco_autocracy2013$location)
mod3.2=coeftest(mod3.2, vcovHAC, cluster= navco_autocracy2013$location)
def3=coeftest(def3, vcovHAC, cluster= navco_autocracy2013$location)
def3.2=coeftest(def3.2, vcovHAC, cluster= navco_autocracy2013$location)
mob3=coeftest(mob3, vcovHAC, cluster= navco_autocracy2013$location)
mob3.2=coeftest(mob3.2, vcovHAC, cluster= navco_autocracy2013$location)


stargazer(mod3,mod3.2, def3, def3.2, mob3,mob3.2, omit =c("location","year"), type="text",
          #title = "Surviving Challenges to Authoritarian Regimes",
          column.labels= c("Challenge Success","Elite Defections", "Mobilization for Autocracy"),
          column.separate = c(2,2,2),
          covariate.labels = c("Minority excluding majority", "Log GDP", "Log population", 
                               "Campaign size", "Foreign support-regime ","Foreign support-campaign",
                               "Opposition unity", "Nonviolent challenge","challenge duration (Log)"
          ),
          no.space = TRUE,
          notes = c("standard errors clustered by country"),
          model.names = FALSE,
          dep.var.caption = "",
          add.lines=list(c("Observations", "1,439 ", "1,188", "1,377", "1,145 ", "1,419 ", "1,168"), 
                         c("Country and Year FE", "No", "Yes", "No", "Yes", "No", "Yes")),
          dep.var.labels.include = FALSE)



# Predicted probability 
mod3= glm(success~ minority, data=navco_autocracy2013, family="binomial")
pp=predictions(mod3, newdata = datagrid(minority = c(0, 1)),vcov = TRUE )
round(pp[2,2]*100, digits = 1) # 0.9
round(pp[1,2]*100, digits = 1) # 6.5

def3= glm(state_defect~ minority, data=navco_autocracy2013, family="binomial") 
pp=predictions(def3, newdata = datagrid(minority = c(0, 1)),vcov = TRUE )
round(pp[2,2]*100, digits = 1) # 3.6
round(pp[1,2]*100, digits = 1) # 11.5




###########################
# Synthetic Control Method 
###########################
#Load data 
dat <- readRDS("dat.RDS")


# Minority excluding a single majority
out_minority <- gsynth(v2caautmob_osp ~ minority+ log_e_gdppc+ log_e_pop+ v2exl_legitideolcr_1+
                         v2exl_legitideolcr_4+ v2x_clphy,
                       data= dat, index = c("ccode","Year"), 
                       force = "unit", CV = TRUE, r =c(1, 5), se = TRUE, inference = "parametric",
                       nboots = 1000, na.rm = TRUE, min.T0 = 12, 
                       parallel = TRUE, estimator = "ife", EM = TRUE)

# plot1
p1 <- plot(out_minority, type = "counterfactual", xlim = c(-12,15),  theme.bw = TRUE)
data_plot <- p1$dat


data_plot$label <- ifelse(data_plot$type == "tr", "Treatment Average", "Counterfactual Average")

plotc1 <- ggplot(data = data_plot, aes(x = time, y = outcome, linetype = label)) +
  geom_line(color = "black") +  # Both lines will be black
  ggtitle("(a) Minority excluding a single majority") +
  xlab("Years since regime onset") +
  ylab("Mobilization for Autocracy") +
  geom_vline(xintercept = 0, linetype = "dotted", color = "#f03b20", size = 0.3) +
  geom_hline(yintercept = 0, linetype = 1, size = 0.2) +
  scale_linetype_manual(values = c( "dashed", "solid")) +  # Make one line solid and one dashed
  theme_minimal() +
  theme(
    plot.title = element_text(size = 13, margin = margin(b = 14), family = "serif"),
    axis.text = element_text(size = 9),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10),
    plot.margin = unit(c(1, 1, 1, 1), "cm"),
    legend.position = "bottom",  # Place legend underneath the plot
    legend.box = "horizontal",  # Display legend items horizontally
    legend.title = element_blank()  # Remove legend title
  )




# Minority excluding other minorities
out_frac <- gsynth(v2caautmob_osp ~ frac_minority+ log_e_gdppc+ log_e_pop+ v2exl_legitideolcr_1+
                     v2exl_legitideolcr_4+ v2x_clphy,
                   data= dat, index = c("ccode","Year"), 
                   force = "unit",CV = TRUE, r =c(1, 5), se = TRUE, inference = "parametric",
                   nboots = 1000, na.rm = TRUE, min.T0 = 12, 
                   parallel = TRUE, estimator = "ife", EM = TRUE)




p2 <- plot(out_frac, type = "counterfactual", xlim = c(-12,15),  theme.bw = TRUE)
data2_plot <- p2$dat

# plot
data2_plot$label <- ifelse(data2_plot$type == "tr", "Treatment Average", "Counterfactual Average")

plotc2=ggplot(data = data2_plot, aes(x = time, y = outcome, linetype = label)) +
  geom_line(color = "black") +
  ggtitle("(b) Minority excluding other minorities") +
  xlab("Years since regime onset") +
  ylab("Mobilization for Autocracy") +
  geom_vline(xintercept = 0, linetype = "dotted", color = "red", size = 0.3) +
  geom_hline(yintercept = 0, linetype = 1, size = 0.2) +
  scale_linetype_manual(values = c( "dashed", "solid")) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 13, margin = margin(b = 14), family = "serif"),
    axis.text = element_text(size = 9),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10),
    plot.margin = unit(c(1, 1, 1, 1), "cm"),
    legend.position = "bottom",  # Place legend underneath the plot
    legend.box = "horizontal",  # Display legend items horizontally
    legend.title = element_blank()  # Remove legend title
  )



#FIGURE 3. Effect of Minority Regime Onset on Mobilization for Autocracy
library(patchwork)
layout=( plotc1 |plotc2 )

pdf("fig3new.pdf",width=8,height=5)
print(layout)
dev.off()


#Table 4. ATTs of Minority Regimes on Mobilization for Autocracy (GSC Models)
out_minority$est.avg # ATTs average
dim(out_minority$wgt.implied) # number of controlled  and treated

out_frac$est.avg # ATTs average
dim(out_frac$wgt.implied) # number of controlled  and treated



###########################################################################################################
###########################################################################################################


